Occupation number and fluctuations in the finite-temperature Bose-Hubbard model 
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We study the occupation numbers and number fluctuations of ultra-cold atoms in deep optical 
lattices for finite temperatures within the Bose-Hubbard model. Simple analytical expressions for 
the mean occupation number and number fluctuations are obtained in the weak-hopping regime 
using an interpolation between results from different perturbation approaches in the Mott-insulator 
and superfluid phases. These analytical results are compared to exact one dimensional numeri- 
cal calculations using a finite temperature variant of the Density-Matrix Renormalisation Group 
(DMRG) method and found to have a high degree of accuracy. We also find very good agreement in 
the crossover "thermal" region. With the present approach the magnitude of number fluctuations 
under realistic experimental conditions can be estimated and the properties of the finite temperature 
phase diagram can be studied. 
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The Bose-Hubbard (BH) model, first studied in the 
late 1980's pj, has long been considered a rather aca- 
demic testing ground for analytical as well as numerical 
approaches in quantum statistical mechanics. Its most 
interesting feature is the existence of a quantum phase 
transition at zero temperature between a superfluid 
and a Mott-insulating phase. With recent advances in 
the experimental techniques of atom trapping and optical 
lattices, the BH model has regained substantial practical 
importance. As first noticed in . the BH Hamiltonian 
applies to bosonic atoms trapped in a deep lattice po- 
tential. Recent experiments with optical lattices Q| have 
spectacularly confirmed the existence of the superfluid- 
insulator transition. Since the Mott-phase is character- 
ized by a well defined occupation number of the poten- 
tial wells, the phase transition has important potential 
applications in quantum information processing JM [|| , or 
Heisenberg-limited matter- wave interferometry both 
of which require optical lattices with regular filling. Al- 
though techniques to eliminate lattice defects have been 
designed @ an important issue for the implementation 
of these applications under experimental conditions are 
fluctuations in the particle number per site at finite T. 
In this work we analyze these fluctuations by studying 
the finite temperature Bose-Hubbard model. 

The zero temperature BH model has been exten- 
sively investigated in the past by methods such as the 
Gutzwiller projection ansatz [jj, the strong-coupling ex- 
pansion EflElElj 113 an d Quantum Monte-Carlo 
[Til ITU Il6t Il7l llS(. as well as with various mean-field 
approaches (see, e.g., PI0,|2(|)- A very powerful nu- 
merical technique in the case of one spatial dimension is 
the density matrix renormalization group (DMRG) plj . 
used to calculate first and second order (i.e., amplitude 
and number-number) correlations [2^J . Less attention 
has been paid to the nonzero temperature properties of 



the BH system. We note here that the quantum phase 
transition exists in the strict sense only at T — 0. For 
finite T the compressibility, i.e. the change of the parti- 
cle number per site with the chemical potential, is always 
nonzero and thus only an approximate Mott phase exists. 
Furthermore the transition between superfluid and insu- 
lating behaviour passes through an intermediate thermal 
region. Recently this thermal crossover has been stud- 
ied within a slave-boson formalism |23j . In this work we 
study the finite-T properties using a perturbative ana- 
lytic approach in the strong-coupling limit, deriving sim- 
ple analytical expressions for the occupation number and 
number fluctuations which cover both quantum and ther- 
mal effects. To verify these results numerically in one 
spatial dimension, we employ a finite-T version of the 
DMRG approach |2l|. We observe good agreement be- 
tween perturbative and DMRG results, including the de- 
scription of the thermal crossover between the Mott and 
insulating phases. 

We consider a d-dimensional cubic lattice of N non- 
linear oscillators characterised by the following Hamilto- 



(1) 



where [i is the chemical potential, M = J^ fc fik is the 



operator of the total number of atoms, and hp. = a\.a,k 



with at, afc being the usual bosonic creation/annihilation 
pairs. The lattice sites are numbered with a d- 
dimensional index fc = {fei,-- - , fed}, while J^ fe means 
summation over all lattice sites. The first term, 



^hop J 



E 

(ki) 



(2) 



on the RHS of is the hopping Hamiltonian, J > 0; 
the sum extends over nearest neighbours and we assume 



2 



cyclic boundary conditions. The second term on the RHS 
of 1|T|1. H.q = H n \ (fife), is a sum of the nonlinear site 
Hamiltonians 



H ni (h k ) = — n k (ri fc 



fin k 

U _ 2 

— (rife - n) + const, (3) 



where n = 
eraging is 

Z~ lr Tt 



[i/U +1/2 is the filling. The quantum av- 
defined in the standard manner as (■ • • ) = 



(...) 



-/3(n-nj<r) 



IS 



where Z = Tr e 

the statistical sum, and [3 = 1/T is the inverse tempera- 
ture. Note that oscillator units with Ti = k-g, — 1 are used 
throughout this paper. 

While elaborate methods are needed if we wish to ex- 
actly pinpoint lobe boundaries or calculate long-range 
correlations, local quantities like the average number of 
atoms and the number fluctuations may be obtained at 
lesser cost. Let us firstly consider the Mott insulator 
phase. To zero order in hopping, the ground state of 
the Mott insulator is |0) J=0 = \\ k |n ) fc , where n — 
round(n), cf. Eq. J2J. A first order perturbative correc- 
tion to this state will contain all states which are found 
from |0)j =0 by moving one atom to a neighboring site. 
All such "single-hop" states are eigenstates of Ho with 
the same relative energy U, resulting in an especially 
simple expression for the ground state to first order in J 
|24|. Up to second order we find 



|0) J=0 + ^(^hop), (4) 



where a is a normalisation factor. By direct calculation 
we then find 



= <0! 



(0 |n*| 0) = no, 

ln\ m\* in\ 2 2zJ 2 n (n 
| 0) - (0 |n fe | 0) = 



(5) 



where z = 2d is the number of nearest neighbours in 
the lattice. These formulae are expected to be a good 
approximation deep in the Mott insulator phase where 
a large energy gap exists and the thermal occupation of 
higher levels can be disregarded. 

It is worth noting that, if staying within the perturba- 
tion approach, these expressions cannot be extended to 
the thermodynamic limit N — > oo. To demonstrate this, 
consider the formula thus obtained for the normalisation 
constant, a: 



l + m (0\Hl p \0) 
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1 

U 2 

{Q\Hl ov \{)) = zJ 2 Nn Q {n + l) 



(6) 



These results hold only if zJ 2 Nn (no + 1) <C U 2 , there- 
fore making the thermodynamic limit impossible. (This 
problem went unnoticed in Ref. |24|.l On the other hand, 
local (on-site) quantities should become independent of 



the size of the system when it becomes large enough 
(thermodynamic limit). In fact, our DMRG simulations 
show that Eq. 10) equally holds if zJ 2 Nn (n + 1) > U 2 . 

In the superfluid phase, one has to account for states 
with the total number of quanta different from Nuq. 
In the limit of zero hopping an arbitrary eigenstate is 



a product state, n) 



n fc in*; 



fe- 



With n close to a 



half-integer, it suffices to keep only two states for each 
site, | no) and \tlq + 1), effectively turning bosons into 
fermions. We shall call the subspace of states with 
iik = no, no + 1, the /-subspace. In this subspace we 
introduce standard fermionic creation and annihilation 
operators, c k c\ + c\c k = S k i. Within the /-subspace, 
(Z,m>0) 



K) a k={c k Vn + l) (n + c k c k ) 



(al)' 4 +m = («o + c{c k ) l (c k Vn + l) 



(7) 



Using these relations one can express any bosonic opera- 
tor in fermionic terms. In particular, 

= j {c k c k + n - nj = AE cj,c fc + const, (8) 

where AE = U (1/2 + no — n). Projection onto the /- 
subspace thus makes the Hamiltonian linear so it can be 
directly diagonalised: (omitting an additive constant) 



E, 



fill, 



(9) 



where <pi = 2n(l - 1)/N, E v = AE - 2J(n + l)coscp, 
and fi are related to c k by a d-dimensional discrete 
Fourier transform. The average on-site particle number 
and number fluctuations are then expressed by the aver- 
age number of "fermions" per site, An = (c k Ck), as 

(rife) = n + An, 

(rife) = ^(«o + 4cfc) 2 ^ = ng + (2n + l)An, (10) 
Sn 2 = (rife) — (rife) 2 = An(l — An). 
For An we readily find: (f3 = 1/T) 



(11) 




This formula holds under the thermodynamic limit N — > 
oo. 

The relations (J5J and l|10|) . (|ll|l are found under mu- 
tually exclusive conditions. Equation J^J follows if we 
assume that only the ground state matters, and retain 
the first non-vanishing correction to its wave function; 
with only one level being important temperature is of no 
concern. Conversely, equations i|l(J|) . (|llf) ignore correc- 
tions to the eigenstate wavefunctions yet account for en- 
ergy shifts; with many levels being accounted for, thermal 
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properties are retained. The question we now address is 
whether a relation can be derived which covers both insu- 
lator and superfluid regions, plus, even more importantly, 
the crossover region (which is naturally termed the ther- 
mal region)? It would be natural to include perturbative 
corrections to all "fcrmionic" states, but this results in 
some intractable algebra. It turns out, however, that for 
J <C U , simple interpolating formulae betwen the two 
above results may be found which apply not only to the 
insulator and superfluid regions, but also to the thermal 
crossover region between them. Namely, 



(n k ) = n + An, 
5n 2 = An(l - An) 



2zJ 2 (n Q + An)(n + An + 1) 



IP 



(12) 



where An is given by Eq. i|llfl . The idea behind l|12|) is 
relatively simple: for small J, the quantum contribution 
(JSJ) is negligible in the thermal and (even more so) the 
superfluid regions, while in the insulator region, An rs 
(or 1), so that (|12l) coincide with JSJ. 

To verify these results in one spatial dimension, we 
use a finite-temperature version of the DMRG approach 
|2l) . At the core of our approach is the following block- 
doubling algorithm: assuming that we know the opti- 
mised states for a block of length L, we can then find 
eigenstates of a ring of three such blocks, (recall that 
we use periodic boundary conditions) and calculate the 
thermal p-matrix. Tracing over the states of one block 
yields the p-matrix of a double-sized block. Diagonalis- 
ing this p-matrix and taking a certain number of eigen- 
states corresponding to the largest eigenvalues results in 
the optimised basis in a block of length 2L. The loss of 
probability is assessed by summing the neglected eigen- 
values; if this loss becomes unacceptable the iterations 
are stopped. To initiate the algorithm, we start from a 
ring of Lq + L' oscillators, and then trace out the states 
of L' oscillators, resulting in an optimised basis for a 
block of Lq sites. We found this to yield much better re- 
sults than starting from an open-ended block by simply 
selecting its lowest eigenstates. The full details of the 
algorithm will be discussed elsewhere. 

A word of caution is necessary here. While borrowing 
much of the technical side of the DMRG method, our 
approach is conceptually different, the difference resid- 
ing in the importance of temperature. Similar to other 
renormalisation-group methods, DMRG relies on the fact 
that the dimension of the physically relevant subspace of 
the full Hilbert space does not grow with the system size. 
For finite temperatures this can only be true until the 



thermal correlation length is reached. Consequently we 
only grow the block up to a size comparable to the ther- 
mal correlation length, and must simply stop as soon as 
the probability loss becomes unacceptable. We use pe- 
riodic rather than open boundary conditions based on 
the same argument. For the temperatures considered 
here the maximum block length reached in this way is 
long enough such that finite-size effects are unimportant. 
For higher temperatures a combination of stochastic and 
DMRG techniques can be used, which we are developing 
and will be discussed in a later work. 

The results of our calculations using both DMRG and 
perturbative approaches in 1-D are summarised in the 
figure, where we plot average on-site numbers and num- 
ber fluctuations for three "cross-sections" along the \i- 
axis, for J = 0.01Z7, 0.02J7, 0.05U. Each plot shows 
data for two temperatures, T = 0.Q1Z7 and T = 0.001C7. 
The top and bottom rows of plots represent, respectively, 
the on-site population, {hu), and number fluctuations, 
Sn 2 = (n|) — (nfc) . To gain more insight into the ther- 
mal region, the middle row of plots shows the difference, 
An, between the on-site population and the nearest inte- 
ger value. At T = 0.001 U, the insulator region is clearly 
defined by An abruptly falling to zero. At T = 0.01U, 
An is a smooth function vanishing as /z goes deeper into 
the insulator region. A similar effect is seen for the on-site 
number fluctuations. At T = 0.001C/, the phase transi- 
tion points are well defined, whereas at T — 0.01J7 the 
boundaries of the insulator phase are "eroded" so that 
Sn 2 is a visually smooth function of u. 

Comparing the DMRG to perturbative results, we see 
that the latter provide a surprisingly good description of 
the quantities in question. For J = 0.01U (left column 
of plots), we find a very good agreement between the 
perturbative and DMRG results. This agreement is rea- 
sonable for J = 0.02C7, (center column) becoming only 
fair at J = 0.05U (right column). Note that, even in the 
latter case, the error is mostly in positioning the ther- 
mal region between the insulator and superfluid phases. 
Away from the thermal region, the results of the pertur- 
bative approach show good agreement with the DMRG 
results, both for the superfluid and for the insulator. 
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